Fracture surfaces of heterogeneous materials: a 2D solvable model 
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Using an elastostatic description of crack growth based on the Griffith criterion and the principle 
of local symmetry, we present a stochastic model describing the propagation of a crack tip in a 2D 
heterogeneous brittle material. The model ensures the stability of straight cracks and allows for the 
study of the roughening of fracture surfaces. When neglecting the effect of the non singular stress, 
the problem becomes exactly solvable and yields analytic predictions for the power spectrum of the 
paths. This result suggests an alternative to the conventional power law analysis often used in the 
analysis of experimental data. 
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The field of fractography is concerned with the char- 
acterization and study of surfaces produced when solids 
break 0. In contrast to other pattern formation mech- 
anisms governed by a phase transition or material de- 
position, in fracture phenomena the broken surfaces are 
created by irreversible material separation induced by the 
propagation of a crack front. A common characteristic 
of many materials is that these surfaces are not smooth. 
Concerning this feature, the pioneering work of Man- 
delbrot el al. and further developments 0, 0, 0, @ 
have shown that crack surfaces can also exhibit frac- 
tal characteristics. For a one dimensional curve, say 
h(x), embedded in 2D geometry, this self-affmc prop- 
erty is characterized by an exponent C, called the rough- 
ness exponent. There are many ways to quantify rough- 
ness Q , one of them being based on the power-spectrum 
P(q) = (h q h_ q ). For a self-affine path, it is expected that 
P(q) ~ g~ 1-2 ^ for intermediate, yet small, q's. Measure- 
ments of the roughness exponent of crack paths in quasi- 
two-dimensional materials have been performed on Berea 
sandstone (C ^ 0.8) 0, concrete (C ^ 0.75) and paper 
(£ ~ 0.6-0.7) @. These variations question the concept 
of universality of fracture roughness exponent in 2D ge- 
ometries. 

A theoretical description of the fracture surface must 
follow from the study of crack tip propagation. A key 
ingredient for that comes from the divergence of the elas- 
tostatic stress field in the vicinity of the crack tip 

^■M)= £ -^^f(9)+T5 ix S jx +0(Vr) , (1) 
£=1,2 y 27Tr 

where (r, 9) are polar coordinates with r = located at 
the crack tip, the x-direction is parallel to the extension 
of the crack path at r = 0, and E^. are known functions 
describing the angular variations of the stress field com- 
ponents [8j. In this expansion, Kg {I — 1,2) and T are 
the static stress intensity factors (SIFs) and the nonsin- 
gular T-stress, respectively. In the framework of linear 
elastic fracture mechanics, the crack propagation is gov- 
erned by the Griffith energy criterion and the Principle 



of Local Symmetry (PLS) |8j. In its simplest form, the 
Griffith criterion states that the crack must grow by sat- 
isfying K\ — Kic, where K\ c is the material toughness. 
The PLS, formulated as K2 — 0, imposes a pure opening 
mode for the stress field in the vicinity of the crack tip. 
Within these two criteria, the crack path is mainly se- 
lected by the PLS, while Griffith's criterion controls the 
intensity of the loading necessary for the crack to grow. 
In addition, the nonsingular T-stress in Eq. rules the 
stability of crack propagation. Stability analysis around a 
straight crack path under pure tensile loading has shown 
that if T > straight crack propagation is unstable, while 
stable when T < Q. 

In this work, we model crack propagation in 2D brittle 
heterogeneous media by using the above elastostatic ap- 
proach formulated in terms of the SIFs and the T-stress. 
We suppose that the crack propagates smoothly by sat- 
isfying the Griffith criterion and the PLS until it meets 
a heterogeneity which induces a change in the fracture 
toughness and/or a local shear perturbation. In response, 
the crack changes locally its direction of propagation by 
forming a kink. If the distance between kinking events 
is small one can show that at leading order the crack ex- 
tensions are straight [10], As shown in Fig. fj] the local 
kinking angle at the position i is given by 59i = 0i+i — 0i, 
which in the continuum limit is S9(x) = h"(x)5x. Now, 
in order to write an equation for the crack path h{x) one 
should relate 89 (x) to the material heterogeneities and to 
the elastic quantities. 




FIG. 1: Modeling crack propagation in a heterogeneous 
medium as a succession of smooth parts and kinks due to 
local defects. 
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Consider a crack in a material subjected to a load- 
ing such that in the absence of defects, a straight line 
h(x) — is the solution for the crack path. For this con- 
figuration, the variation of the kinking angle 89 is related 
to the variation of the mode II stress intensity factor 8K2 
induced by the defects 0, 0] . To linear order, this re- 
lation is given by 86 = -2{8K 2 )/K X [nl[ll|. Next, one 
can write SK2({h},x) as a functional of h(t), the shape 
of the crack path for t < x 12] . Putting all together, 
one obtains the following dimensionless equation for the 
crack path 



h"(x)=r l2 (x)+r ] [(x) 



h'{x) + ah(x) - (3 



h'{t)dt 



\Jx — t 



(2) 

A systematic derivation of this equation will be given in 
[l3| . In Eq. @, the primes denote derivatives, a is a 
positive quantity fixed by t he g eometry of the sample in 
which the crack propagates [12j and f3 is proportional to 
the T-stress. In addition, Eq. J2J) contains two space de- 
pendent functions which represent the local fluctuations 
in the toughness, modeled by 771(2;), and local shear fluc- 
tuations, modeled by 772(2;)- Note that when the material 
is homogeneous, f]' e (x) =0 (£ — 1,2), and the solution of 
Eq. J2Jl is given by a straight path, recovering the zeroth 
order solution. Therefore Eq. Q describes the devia- 
tion from a straight crack propagation in the presence of 
inhomogeneities. 

The numerical resolution of Eq. (0) reveals that when 
> the crack paths are unstable, while they are stable 
for f3 < ^3 ■ This is consistent with the T-criterion || , 
however we found that for p < the crack propagation 
is more stable in our formulation, since the crack path 
decays to a straight line rather than exhibits a square- 
root behavior (h(x) ~ y/x) [9j- In the following, we focus 
on the case (3 = for two reasons. First, power counting 
shows that as long as the path is stable the /3-term is 
subdominant at large scales. This is verified by numerical 
simulations and turns out to be valid at all scales [l3| . 
Second, the model is exactly solvable for this case, which 
allows for conclusive predictions. The discrete version of 
Eq. (2J with P = yields 

h i+ i = 2h t - hi-x + (r]2 : i - T/ 2: i-i)A 

+(Vi,i - Vhi-i) {hi - hi-i + ahiA) , (3) 

where 1 < i < N. The quantities rj^i are modeled 
as independent random numbers, characterized by their 
variances (r\i^r\ m .j) = D^8f m 8ij. The initial conditions 
are chosen such that h{x) = for x < 0, implying 
ho = hi = 0. In addition, the x-axis is scaled such that 
the observation window is always < x < 1, thus for a 
ID lattice with N sites, one has A = 1/N. Fig. |2shows 
an example of a path grown using this procedure and 
shows that the model produces rough crack lines with 
possible self-similar properties. 




FIG. 2: An example of a crack path produced using Eq. 0, 
with A = 10~ 5 , a = 20, Di = 0.04 and D 2 = 6.25. The inset 
shows a zoom of a smaller part of the path. 



Let us first study crack paths when there is a single 
shear perturbation given by 772. i = #o<5i2- This situation 
is equivalent to redefining the boundary conditions as 
h(0) = and h'(0) = 80 [9|- Later, we will generalize our 
results to any uncorrelated shear perturbations 772,1 and 
find that this source of randomness is essential to repro- 
duce the experimentally observed rough crack paths. For 
V2,i = 0{)8i2i we are going to show that crack paths tend 
to decay into straight ones (see insets in Fig. 01) • This 
decay can be also accompanied by oscillations. Interest- 
ingly, when considering an ensemble of paths grown using 
the same parameters 8q, a and Di/A, but with different 
choices of the noise 771 one finds that for small A, the 
solution becomes independent of the realization of the 
noise. To understand this, let us consider the logarith- 
mic derivative of h, namely ipi — hl ^ h h ^ 1 14]. In terms 
of ipi, Eq. © with r)2,i = can be rewritten as 



i+l 



a (77; 



i.i - 77i,j-i) + 



j>i + (7?i,i - 77i,i_i) ipj 
1 + ViA 



(4) 



From the observation that the solutions depend on D\ / A 
only, 771 should be of order 0{\fK). We will be interested 
in deriving a Fokker-Planck equation for the t/>'s. As 
usual in stochastic equations, the analysis of Eq. for 
small A requires keeping terms up to order A, since the 
leading order (order \/A) has zero average. This gives 



i+i 



4>i + im.i - vi,i-i) + a ) ~ ^ 2 A ■ (5) 



Next, using Eq. (Jjjjl one writes the integral equation for 
the probability distribution of the t/>'s, namely P(^,r]i), 
which is readily given by 



Pi+i (ipi+i,t]i) = J d77, i _ 1 p(?7 i _ 1 ) / dipiPi{ipi,rn-i) 
8 (ipi+i -ipi- fa - 7/j_i) (ipi + a) + A) , (6) 

where p(r)i) is the distribution of the random toughness 
at the location i. For brevity, the subscript 1 in 771 has 
been suppressed. In order to integrate over ipi and to 
eliminate the 8 function from the integral in Eq. 10, one 
needs to express ipi as a function of tpi+i . To achieve that 
we iterate Eq. (J5J to first order in A and then get 
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ypi+i - (m - w-i) (i>i+i + (vi - m-i) 2 (^+1 + a ) + ^-1) 



Pi+i (ip i+ i,Vi) = / dr)i^ip(r)i-i) — — r — — . (7) 



I 

The denominator in Eq. (J2J) results from the Jacobian of the (^-function. Expanding the integrand to order A gives 
I 



P i+ i(tp,rn) = J drii^ipirii-!) ^Pity,^-!) - fa - rii-i) — [(4> + a) Pi(i),r)i-i)] 
+ (m - Vi-if ^ [l (V> + «) 2 Pi Vi-i)] + A^; [^P* (V, Vi-i)] 



(8) 
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We see that, at order A, the probability distribution 
function Pj+i is quadratic in 77.;. Thus one can write 

Pii^m-i) = P^W+vt-iP^W+vliP^W • (9) 

Plugging this expansion into Eq. JSJl and 
performing the integral over 77.^1 gives 



d 

' dtp 



ty + a)if>(^) and P£[ (0) 



>(2) 



^2 | (V» + a) 2 P^ ("0) ■ an d finally we get from 
Eq. (jHJ the x-dependent Fokker-Planck equation for p(°> 



dx 



dip 



(i> 2 + C(tP + a)) P(°) 



(10) 



where C = with D x = (r)f). Note that Eq. HH does 
not have second derivatives of P^ with respect to ?/>. 
Thus, strictly speaking, it looks like a Liouville equation. 
This means that the probability distribution function of 
ip becomes that of a deterministic evolution given by 



P i+1 (tp i+ l) = J dipiPiiipt) 

S (^+i -1>i+ + Ctyi + a)] A) , (11) 

from which we can read the deterministic equation 

d _± = _ (/ . a 
dx 



-4> 2 -C(ip + a) . 



(12) 



Since in the continuum limit ip(x) = h'(x)/h(x), Eq. i|12|) 
yields 



h"(x) = -C[ti{x) + ah(x)\ 



(13) 



Thus, one finds that the effective equation (|13|) is ob- 
tained from the noisy one (Eq. J2J with = 0) by sim- 
ply replacing the noise term r)[(x) by a negative con- 
stant, — C, proportional to its variance, even though 



jfi(x) is equally positive and negative. This nontrivial 
result shows that, in the small A limit, the local tough- 
ness fluctuations r/i^ do not influence the shape of the 
crack path, apart from setting a relevant scale for the 
energy that has to be invested in making the crack prop- 
agate. At this point, one can easily solve Eq. (|13fl with 
the initial conditions h(0) = and h'(0) = 9q: 



h(x) = 9$e 



Jx sinh (±VC 2 -AaCx) 



AaC 



(14) 



The solution for the crack path (jl4|l exhibits either an ex- 
ponential decay or damped oscillations depending of the 
sign of (C — 4a). For the present case, the power spec- 
trum of the crack path can be obtained analytically. It is 

r 2 i _1 

given by (h q h_ q ) = 6 2 (q 2 — Co) + C 2 q 2 , where h q 

is the Fourier component of h(x). In Fig. |3J we compare 
this expression for (h q h- q ) with the corresponding aver- 
aged power spectrum over ten realizations of the noise 
771. As can be seen, both the cases C < 4a (damped 
oscillations) and C > 4a (exponential decay) show that 
the theoretical results fit the numerical data. 

Our next step is to consider crack propagation when 
there are many shear perturbations. This amounts to 
retaining the additive noise rj 2 in Eq. ©. Unlike the 
fluctuations in the local toughness 771, the shear pertur- 
bations cannot be modeled by a deterministic term. Fur- 
thermore, its randomness seems crucial for creating the 
rouehpatterns similar to those observed in experiments 
[IS E| (see for example Fig.EJ. Interestingly, by vary- 
ing the amplitude of the noise D2, one can can produce 
rather different patterns. Based on the previous results 
we generalize the expression for the power-spectrum to 
the current case. Note that in Eq. ©, 772 is just a non 
homogeneous term. Thus, once we have a solution for the 
homogeneous equation, we can construct a special solu- 
tion for the non homogeneous one. Using the fact that 771 
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teresting to use expressions like Eq. (1151) instead of sim- 



L °3',o« 



Log, q 



FIG. 3: Power spectra of crack paths with a single shear 
perturbation averaged over 10 realizations, using 9o — 1 and 
A = 4 x 1(T 4 . (a) a = 100, Di = 4 x 10" 3 (C < 4a). (b) 
a = 20, Z?i = 4 x 10~ 2 (C > 4a). The solid lines correspond 
to the theoretical curves. The insets show examples of such 
paths. 




FIG. 4: Power spectrum of crack paths with many shear per- 
turbations averaged over 10 realizations, using A = 10~ 4 , 
a = 1, D\ = 0.04 and D2 = 4.2. The solid line corresponds 
to the theoretical prediction. The inset shows an example of 
such a path. 



and 772 are independent random variables, we conclude 
that averaging over realizations of the local toughness 
fluctuations still amounts to replacing Tj[(x) by — C. Do- 
ing so, we write the power spectrum for a general noise 
7/2(2;) in the form 

(l + C'-^q 2 + (2C + a + 2) 

w -»VU« ' (15) 

which is different from the power-spectrum in the case 
of a single shear perturbation, as there is a q 2 term in 
the numerator. The coefficients in Eq. (|15|l are deter- 
mined from the Fourier transform of a stationary sig- 
nal after cutting out the transient regime, that is h q = 
J Xg h(x)e lqx dx. This leads to a system with effective ran- 
dom initial conditions at a; = xq. Since xq is chosen in 
the steady regime the statistics of h(xo) and h'(xo) are 
known, and an analytical expression of the power spec- 
trum can be obtained. Fig. 0] shows that despite Eq. i|TB|) 
agrees very well with the simulations data for the power 
spectrum, one might be tempted to fit it with a power 
law ansatz. However, apart from a tail q~ 2 at large q's, 
which yields a roughness exponent £ = 0.5 at small length 
scales, Eq. 1)150 teds us that there is no self-affine behav- 
ior of fracture surfaces at intermediate length scales. 

Here we have introduced a model that describes the 
motion of a crack tip in a 2D heterogeneous brittle mate- 
rial and allows for the analysis of the form of the surface 
created. The approach is based on a pure elastic de- 
scription of the crack growth and does not necessitate 
additional ingredients such as propagation through voids 
nucleation 15]. The results for the stability of a straight 
crack propagation are in agreement with the well known 
T-criterion [9j. Exact results are provided for the case 
(3 = 0, that allow to predict analytically the power spec- 
trum of the crack surface. These results suggest that the 
roughness of the surface is not self-affine, but rather a 
fuzzy shape convoluted with an oscillating path. 

From an experimental point of view, it could be in- 



ple power-laws to fit the data of rough crack lines in 
2D geometries. However, such analysis require averag- 
ing over many repetitions under the same experimental 
situations, which is difficult to perform due to the irre- 
versibility of fracture events. Finally, our conclusions are 
not relevant to crack surfaces only, but concern analysis 
of rough shapes where an oscillating background is likely 
to exist and might bias the measurement of the surfaces 
characteristics |7J, LU| • 
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